rm(list = ls())

library(foreign)
library(MASS)
library(AER)
library(splines)
library(mgcv)
library(boot)
library(rgl)
library(margins)
library(lme4)
library(survival)
library(coxme)
library(countrycode)
library(readstata13)


setwd("/Users/davidcarter/Dropbox/Systemic Crisis Rep/")


Data<-read.dta("data_w_allmeasures_int.dta")


Data<-na.exclude(Data)
attach(Data)


#############################
#
####### Regular PCA #########
#
#############################


fit_hg1<-gam(make_claim_onset_f ~number_of_changes+ pca_1 +ti(number_of_changes, pca_1)+factor(nosettle)+latitude+longitude+shape_area_10+s(onsetyrs,k=10),family="gaussian")


number_of_changes1= matrix(seq(min(number_of_changes,na.rm=T),max(number_of_changes,na.rm=T),length.out=20),20)
number_of_changes2=as.matrix(rep(number_of_changes1,20))
number_of_changes3= matrix(number_of_changes2,length(number_of_changes2),1)
number_of_changes4=matrix(number_of_changes2, 20, 20)


pca_1_1 = matrix(seq(min(pca_1,na.rm=T),max(pca_1,na.rm=T),length.out= 20), 20,1)
pca_1_2= as.matrix(rep(pca_1_1, 20))
pca_1_2 = as.matrix(sort(pca_1_2))
pca_1_3=matrix(pca_1_2,length(pca_1_2),1)

pca_1_4=matrix(pca_1_2, 20, 20)



newd_hv1<-data.frame(cbind("number_of_changes"= number_of_changes3,"pca_1"= pca_1_3,"number_of_changes"= number_of_changes3,"pca_1"= pca_1_3,"nosettle"=0,"latitude"=mean(latitude,na.rm=T),"longitude"=mean(longitude,na.rm=T),"shape_area_10"=mean(shape_area_10,na.rm=T),"onsetyrs"=mean(onsetyrs,na.rm=T)))



colnames(newd_hv1)<-c("number_of_changes","pca_1","number_of_changes","pca_1","nosettle","latitude","longitude","shape_area_10","onsetyrs")



pred3d_hv<-predict.gam(fit_hg1,type="response", newd_hv1,se.fit=T)

pred3d_hv$fit[pred3d_hv$fit <0] <- 0


P3d_hv<-matrix(pred3d_hv$fit,20,20)


persp3d(pca_1_4, number_of_changes4 , P3d_hv,col="gray",xlab = "Systemic Instability", ylab = "Historical Border Variability", zlab = "Pr(Claim Onset)")


## Produce Figure 5a in the Manuscript ##

persp(number_of_changes1, pca_1_1, P3d_hv,col="white",ylim=range(pca_1_4),xlim=range(number_of_changes4),zlim=c(0,0.65),ylab="\nComposite Systemic Instability",xlab="\nHistorical Border Precedents",zlab="\nPr(Claim Onset)",theta=325,phi=11,scale=T,ticktype="detailed",shade=.5)



#################################################
#									   			#
#  Explore Weighted PCA measure of instability  #
#									  			#
#################################################


fit_hg1<-gam(make_claim_onset_f ~number_of_changes+dw_pca +ti(number_of_changes, dw_pca)+factor(nosettle)+latitude+longitude+shape_area_10+s(onsetyrs),family="gaussian")



number_of_changes1= matrix(seq(min(number_of_changes,na.rm=T),max(number_of_changes,na.rm=T),length.out=20),20)
number_of_changes2=as.matrix(rep(number_of_changes1,20))
number_of_changes3= matrix(number_of_changes2,length(number_of_changes2),1)
number_of_changes4=matrix(number_of_changes2, 20, 20)

dw_pca_1 = matrix(seq(min(dw_pca,na.rm=T),max(dw_pca,na.rm=T),length.out= 20), 20,1)
dw_pca_2= as.matrix(rep(dw_pca_1, 20))
dw_pca_2 = as.matrix(sort(dw_pca_2))
dw_pca_3=matrix(dw_pca_2,length(dw_pca_2),1)

dw_pca_4=matrix(dw_pca_2, 20, 20)

newd_hv1<-data.frame(cbind("number_of_changes"= number_of_changes3,"dw_pca"= dw_pca_3,"number_of_changes"= number_of_changes3,"dw_pca"= dw_pca_3,"nosettle"=0,"latitude"=mean(latitude,na.rm=T),"longitude"=mean(longitude,na.rm=T),"shape_area_10"=mean(shape_area_10,na.rm=T),"onsetyrs"=mean(onsetyrs,na.rm=T)))

#newd_hv2<-data.frame(cbind("number_of_changes"= number_of_changes3,"dw_pca"= dw_pca_3,"number_of_changes"= number_of_changes3,"dw_pca"= dw_pca_3,"dyad_m"=19,"nosettle"=2,"onsetyrs"=mean(onsetyrs,na.rm=T)))

colnames(newd_hv1)<-c("number_of_changes","dw_pca","number_of_changes","dw_pca","nosettle","latitude","longitude","shape_area_10","onsetyrs")

#colnames(newd_hv2)<-c("number_of_changes","dw_pca","number_of_changes","dw_pca","dyad_m","nosettle","onsetyrs")

#pred3d_hv<-predict.gam(fit_hg1,type="response", newd_hv1,se.fit=T)

pred3d_hv<-predict.gam(fit_hg1,type="response", newd_hv1,se.fit=T)

pred3d_hv$fit[pred3d_hv$fit <0] <- 0

P3d_hv<-matrix(pred3d_hv$fit,20,20)


persp3d(dw_pca_4, number_of_changes4 , P3d_hv,col="gray",xlab = "Systemic Instability", ylab = "Historical Border Variability", zlab = "Pr(Claim Onset)")


## Produce Figure 5b in the Manuscript ##

persp(number_of_changes1, dw_pca_1, P3d_hv,col="white",ylim=range(dw_pca_4),xlim=range(number_of_changes4),zlim=c(0,.65),ylab="\nComposite Regional Instability",xlab="\nHistorical Border Precedents",zlab="\nPr(Claim Onset)",theta=325,phi=11,scale=T,ticktype="detailed",shade=.5)


